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Abstract. We study a reaction-diffusion equation with an integral term describing 
nonlocal consumption of resources. We show that a homogeneous equilibrium can lose its 
stability resulting in appearance of stationary spatial structures. It is a new mechanism of 
pattern formation in population dynamics that can explain emergence of biological species 
due to intra-specific competition and random mutations. Travelling waves connecting an 
unstable homogeneous equilibrium and a periodic in space stationary solution are studied 
numerically. 



1. Introduction 

The concept of degeneracy in biology was proposed by Edelman in 1978 in relation with 
natural selection in biological systems (see [1]). Degeneracy, according to his definition, 
means the ability of elements that are structurally different to perform the same function or 
yield the same output. He points out that degeneracy is necessary for natural selection where 
natural selection is understood here not only in the context of the evolution of biological 
species but for many other biological systems including cell populations, immune systems, 
brain cortex and others. Degeneracy is related to robustness of biological systems and to 
their ability to compensate lost or damaged functions. For example, some proteins thought 
to be indispensable for human organisms can be completely absent in some individuals (see 
[1]). There are also many examples related to the cortex activity where some functions of 
the damaged area can be taken by some other areas. 

In [2] , Kupiec and Sonigo emphasize many physiological situations such as cell differenti- 
ation or immunology, where darwinian selection processes arise. They show that degeneracy 
and competition are very common phenomena leading to the emergence of structures. 

Atamas studies the degeneracy in terms of signals and recognizers [3]. Each recognizer can 
receive not only the signal, which corresponds it exactly, but also close signals. He suggests 
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that each element recognizes a corresponding signal with a very high specificity, which is not 
however absolute. The same element can recognize other more or less similar signals though 
with less probability. For example, it can be antigen receptors and antibodies that have 
relative but not absolute specificity of recognition. He fulfils numerical simulations with 
the individual based approach where the recognizers can receive signals with degenerate 
specificity and can reproduce themselves. This can be considered as a model for natural 
selection with degenerate recognition. The results of the numerical simulations show that 
initially localized population of recognizers can split into two subpopulations. Therefore the 
degeneracy can result in the emergence of new species due to intraspecific competition and 
selection. 

In this work we study degeneracy in biological system with continuous models. We con- 
sider them in the context of population dynamics. Starting with the classical logistic equation 
describing reproduction and migration of individuals of a biological species, we introduce a 
nonlocal consumption of resources that can influence the death function. One of possible 
interpretation of this model is as follows. Suppose that each individual is characterized by 
his positions {x, y) at time t. The characteristic time scale is related to reproduction and can 
be measured for example in years. The individual can slowly migrate changing from time to 
time his position. For example it can be related to a new nesting place every year or every 
several years. The reproduction takes place at the same space point. However, contrary to 
conventional models in population dynamics we consider the space point {x, y) not as the 
exact position of the individual but as his average position. This means that his nesting 
(reproduction) place is still a point but he can consume resources in some area around the 
point {x.y). Therefore the death function, which takes into account the competition for 
the resources, depends not only on the individuals located at the point {x, y) but also on 
those who are in some area around. We can ask the question, why the individuals coming to 
consume the resources at the point (x, y) influence the death function at this point. One of 
possible answer could be for example that young individuals who are attached to their nest- 
ing place and who cannot yet consume the resources in a bigger area will have less chances 
to survive. We will specify the model in the next section. 

We will discuss also another interpretation of the model related to the evolution of bio- 
logical species. In this case the space variable corresponds to some morphological character- 
istics and diffusion to a random small change of this characteristics due to mutations. The 
degeneracy in this case means that individuals with different values of this morphological 
characteristics can consume the same resources and are in competition for these resources. 
It can be for example the size of the beak for birds eating the same seeds or the height of 
the herbivores eating the same plants. 

We will show that the homogeneous equilibrium stable in the case of the usual logistic 
equation can become unstable in the case of nonlocal consumption of resources. It is a 
new mechanism of pattern formation where spatial structures can appear not because of the 
competition between two species as in the case of Turing structures but for a single species. 
In the context of population dynamics it can be intepreted as emergence of new species 
due to the intra-speciflc competition and random mutations. If the population is initially 
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localized it will split into several separate subpopulations. Its dynamics resembles Darwin's 
schematic representation of the evolution process. Such behavior cannot be described by the 
logistic equation without nonlocal terms. 

Similar to the conventional reaction-diffusion equation, reaction-diffusion equation with 
integral terms can have travelling wave solutions. Depending on the parameters it can be 
either a usual travelling wave with a constant speed and a stable homogeneous equilibrium 
behind the wave or a periodic wave with a periodic in space solution behind the wave. Their 
speeds of propagation (average speed for the periodic wave) are the same as the minimal 
speed for the monostable reaction-diffusion equation. It is determined only by the diffusion 
coefficients and by the derivative of the nonhnearity at the unstable equihbrium. 

We note that there is a number of works where the diffusion term is replaced by an 
integral term describing a nonlocal interaction (see, e.g., [4]). In spite of some similarity of 
the models they are quite different. The scalar equation with the nonlocal diffusion term 
satisfies the comparison principle and cannot describe emergence of spatial patterns and 
propagation of periodic waves. 

Mathematical properties of some integro-differential evolution equations are studied in 
[5,6]. 



2. Model with a nonlocal interaction 



We begin with a well established in population dynamics model taking into account 
migration of individuals. We suppose that they move randomly during their life time and 
all direction of motion are equally possible. If p{r)dr is the probability of the displacement 
on the distance between r and r + dr ai the unit time measured in generations, then 



p = y y r'^p{r)dr 

is a meansquare displacement during one generation (see, e.g. [7]). Assuming that migration 
is independent of the rate of birth and death, we obtain a local change of the density N of 
the population at the point (x, y) during the time interval Ai: 

AiV(a;, y,-t)=(^J N{x', y', t)^dx'dy' - N{x, y, t) + F{N, x, y)^ At, 

where r = {x' — x)"^ + {y' — y)"^. Using the Taylor expansion for N(x',y',t) around {x,y) 
we obtain the nonlinear diffusion equation 

where d = is the diffusion coefficient, F is the local rate of growth of the population. 
Usually it is considered in the form 

F^{B-D)N, (2.2) 



3 



where B is the birth function and D is the death function. The birth function B{N) has a 
sigmoidal form and it is constant for sufficiently large. The death function D is usually 
considered in the form 

D{N) = b + kN, (2.3) 

where the first term in the right-hand side describe the natural death, and the second term 
its increase because of the competition for limited resources. This form of the death function 
implies that the resources are consumed locally, that is exactly at the same point (x, y) where 
the individual is located. 

In a more general and more realistic case the resources can be consumed in some neigh- 
borhood of this point. We take this into account in the death function: 

D{x,y,t) =b + k [ ([ (j){x-x',y-y',t-t')N{x',y',t')dx'dy'^dt'. 

Jo \J{x,y) / 

Here 0(x — x',y — y', t — t') is the function that shows how the individuals located at (x', y') 
and at time t' infiuence the resources at the point (x, y) at time t. This expression for the 
death function means that survival of individuals at the point (x, y) depends on the resources 
available at this point. On the other hand, these resources are consumed not only by the 
individuals that have their average position at this point but also by the individuals from 
some area around it. 

If the motion of the individuals is fast with respect to the given time scale, then the 
location of the individual with an average position at the origin is time independent and 
given by the probability density function (j)Q{x,y). If we assume moreover that the resources 
are renewed with a small characteristic time compared with other time scales (reproduction, 
death, migration), then the infiuence function 0(x — x',y — y', t — t') is concentrated at t = t' 
and equals zero for t' < t: 

(f){x -x',y- y', t-t') = (poix -x',y- y)5{t - t'). 

Thus we obtain the modified expression of the death function: 

D{x,y,t)=b + k ^o{x-x',y-y')N{x',y',t)dx'dy'. (2.4) 

We note that if </)o(x,?/) is the Dirac 5-function, then it is reduced to ()2.3|) . Substituting 
fl2.2|) and ()2.4|) into ()2.ip we obtain the integro-differential equation 

ON fd'^N d'^N\ f f \ 

We will analyze this model in the one-dimensional spatial case where B{n) is a constant. 

3. Evolution of species 

In the previous section we considered a spatial distribution of individuals in a population 
and its evolution in time. Let us characterize the population not by individuals but by genes. 



4 



Consider a gene G and suppose that due to mutations it can be in different states, which we 
denote by Xj, i = 1, A^. Let Ciit) be the number of genes in the population at the state 
Xi at time t. We suppose further that mutations can happen only between the neighboring 
states, that is for the state Xi the gene can go to one of the states Xi-i and Xj+i. If we 
assume that all transitions are reversible and have equal probability p, then the number of 
genes changing their state from to Xj during the unit time will be pcj_i. The number 
of opposite transitions is pcj. Therefore the flux of the concentration at Xi is 

fi = J9(q_i - Ci). 

The probability p of the transition depends on the distance between the states and will be 
considered in the form 

d 

p = — z — ' 

where d is some constant. We suppose for simplicity that the distances between the neigh- 
boring states are the same and that d does nor depend on i. The last expression is well 
defined if the distance between the states is sufficiently large. 

Thus, under the assumptions above we come to the usual diffusion fluxes. Passing from 
the discrete to continuous model we obtain the equation 

Be cP' c 

(3.1) 

where c(x, t) is the concentration of genes at the state x at time t, and, as above, i? is a 
birth function, D is a death function. Taking into account the nonlocal consumption of the 
resources we come to the same model as in Section 2. 



4. Stationary solutions 

Consider the equation 

= J ^^^^ y)(^hj)dy\ - he (4.1) 

for X G M. Here a and h are some constants, is a function with a bounded support and 
such that 

/>oo 

= 1. (4.2) 

We will assume that (/)(?/) is nonnegative, even and identically equal zero for \y\ > N. If its 
support tends to zero, then we obtain in the limit the 5-function, and equation 1)4.11) becomes 
the usual reaction-diffusion equation 

5f = <'a;i + '^('>-c)-fc (4.3) 



often considered in population dynamics and in other applications. 

It can be easily verified that equation ()4.H) has two homogeneous in space stationary 
solutions, co(x) = and ci{x) = a — b. These are also stationary solutions of equation ()4.3|1 . 



5. Stability 

If a > 6, then Ci{x) is a stable stationary solution of equation ()4.3p . We will show that 
it can be unstable as a solution of equation ()4.1|) . Linearizing this equation about Ci, we 
obtain the eigenvalue problem 

POO 

dc — a (f){x — y)c{y)dy = Ac, (5.1) 



where a = a — b. Applying the Fourier transform to this equation, we have 

(rfe' + ^0(O + A)5(O=O, (5.2) 

where the tilde denotes the Fourier transform of the corresponding function. This equation 
has nonzero solutions at the points of the spectrum, that is for such values of A that 

d^'^ + a0(O + A = 0. (5.3) 

Consider the function 

If it becomes negative for some values of C,, then there are positive A satisfying equation 
()2.4|1 . In this case, the stationary solution ci is unstable. 
We note that 

m = / mdy = 1. 



Therefore, for any given continuous and bounded function 0(0; ^{0 becomes strictly posi- 
tive if d is sufficiently large. Thus, there is no instability if the diffusion coefficient is large. 
By virtue of the assumption that the function (f){y) is even, we have 

/oo 
(p{y) cos{^y)dy. 
-co 

The assumption that the function (j){y) has a bounded support and condition ()4.2|1 are not 
essential for the stability analysis. For the functions 

My) = e-''y\ 02(y) = e""l^l (a>0). 



we have 



MO = J- e-^''^'^\ MO 
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respectively. Therefore, $(^) is positive for all ^, and the stationary solution ci is stable. 



6. Piece-wise constant distribution 

Consider the following function: 



Then 



1 1 
0(0 = ^ y ^cos{^y)dy = — sin(eiV) 

m = de + ^ sin(eiV). 

Depending on the parameters, this function can be strictly positive or can change sign. It 
can be easily verified that $(0) = a > 0, and it is positive for |^| sufficiently large. If the 
diffusion coefficient is large, then $(0 is positive for all ^, and the stationary solution is 
stable. The instability can appear if d is small enough. 

In the critical case, in which the function is nonnegative and not strictly positive, 
it satisfies the equalities 

$(eo) = o, $'(eo) = o 

for some ^o- K follows from these relations: 

+ 7^ sin(eoiV) = 0, 2d^o - tSt? M^oN) + f cos(eoiV) = 0. 

Let us introduce the following notations: = d/a, z = ^qN. Then the previous equations 
can be written in the form 

^9 sin 2; r, sin 2; 

/i^o + = 0, 2/ieo + cos z = 0. (6.2) 

z z 

Eliminating /x^q, we obtain from these two equations the equation with respect to z: 

tan 2; = -z. (6.3) 

This equation has an infinite number of solutions. We will consider only positive solutions. 
Denote them by zi,Z2, ■■■ For each of them we can find from the first equation in ()6.2|) the 
relation 

sin Zi 

/i = -iV^, J = 1,2,... (6.4) 

Thus we have a countable number of curves on the {N, /x)-plane. Each of them gives the 
values of parameters for which the spectrum of the eigenvalue problem ()2.H) crosses the 



and 
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origin. We note that sinzj is positive for even j and negative for odd j. Since /i is positive, 
we will consider only odd values of j. 
It can be easily verified that 

sin Zl sin 

z\ zl 

Therefore 

/ii(Ar) > ii-^{N) > ... for N>0. 

We have seen in the previous section that there is no instability if d is sufficiently large. 
Therefore, fi > yUi(A^) corresponds to the stability region, fi < fii{N) corresponds to the 
instability region. The stability boundary is given by the equation 

/i = /ii(iV). (6.5) 

We can also find the period r of the spatial structure at the stability boundary: 

_ 27r _ 27rN 
Co Zl 

Let us consider a fixed value of n and vary A^. For small enough the homogeneous 
stationary solution is stable. For some critical value of such that ()6.5|) is satisfied, it loses 
its stability. For greater values of A^, the solution remains unstable and other points of the 
spectrum can cross the origin. 

Suppose now that we consider a bounded interval with the length L and with the periodic 
boundary conditions. The analysis above remains valid but we need to take into account 
the period of the perturbation. We fix ^ = 27i/L and consider the equation $(^) = 0. We 
obtain the equation sin z = —fi^'^z. If n is sufficiently large, then the equation does not have 
nonzero solutions. In the critical case we should also equate the derivatives of the function 
in the left-hand side and in the right-hand side of the last equality: cos 2; = — yU^^. Therefore 
we can find z from the equation tanz = z, and the critical value of fi from the relation 
/i = -{smz)/{^^z). 



7. Emergence of structures 

Equation ()2.1|) is considered on the segment [0, L] with periodic boundary conditions. 
It is discretised by using an explicit finite differences scheme. The function $ is piece-wise 
constant as described in equation ()6.H) . The integral is approximated by the trapeze method. 

Figure 1 shows the behavior of the solution of equation ()4.H) for the values of parameters 
d = 0.05, a = 2, A^ = 3, 6=1 and L = 10 and for the initial condition c{x, 0) = a — b 
perturbed by a small random noise. This initial homogeneous equilibrium is unstable and 
a spatial structure appears. Note that in accordance with the stability analysis (Section 6) 
the homogeneous equilibrium is unstable only for small values of the diffusion coefficient d. 
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Figure 1: Emergence of spatial structure 




Figure 2: Splitting of a population 
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Figure 3: Evolution 

For the same values of the parameters above and d > 0.11, the homogeneous equilibrium is 
stable. 

We consider next the same values of parameters as for Figure 1 and the initial condition 
c(0, x) = a — 6 in a small interval centered at the middle of [0, L] and c(0, x) = otherwise. At 
the first stage of its evolution the solution grows and diffuse slightly (Figure 2), its maximum 
remains at the middle of interval [0,L]. At the second stage the solution decreases at the 
center and grows near the borders of the interval. We observe here that the intra-specific 
competition results in appearance of two separated sub-populations symmetrically situated 
around the middle of [0, L]. 

The case of a larger interval with two initially separated sub-populations is shown in Fig- 
ure 3. Each of them first splits into three sub-populations. Those in the center do not change 
any more. The sub-populations from the sides split once more into two sub-populations each 
of them. These consecutive splitting into sub-populations has an interpretation related to 
the evolution of species. We discuss this question below. 

The numerical simulations shown in Figure 4 are carried out for the same values of the 
parameters except for the diffusion coefficient which is smaller in this case. The consecutive 
sub-populations are almost disconnected here. Each new structure emerges from very small 
values of the concentration. 

We consider finally the case where the function (j){y) is not symmetric. It is similar to 
that given by ()6.H) but equals for ?/ < (see Section 9). The dynamics of the solution is 
different in this case: the first maximum of the concentration moves to the left, then the 
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Figure 4: Small diffusion coefficient 

second maximum appears and moves remaining at fixed distance from the first one, and so 
on (Figure 5). 



8. Approximate solution 

In this section we consider an approximation of equation ()4.H) . which is not well justified 
mathematically but allows us to explain some of its properties. Consider the integral in the 
right-hand side of this equation 

/oo 
(f){y)c{x - y)dy 
-oo 

(the time dependence of the function c is for brevity omitted). We will use the Taylor 
expansion 

c(x — y) = c{x) — c {x)y + -c {x)y'^ + ... 
If (j) is an even function then we obtain 

0(a; - y)c{y)dy = c{x) + 7c" (x) + (8.1) 



where 



7=2' (t>{y)y'^dy- 
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Figure 5: Asymmetric evolution 



Taking only the first term in the right-hand side of expansion ()8.1|1 . we obtain the reaction- 
diffusion equation ()4.H|1 instead of the integro-differential equation ()4.1|1 . Equation ()4.Hj) is 
well-studied. It is known that it has a travelling wave solution, that is the solution of the 
form c{x, t) = u{x — st), where s is the wave velocity. In fact, monotone in x waves for this 
equation exist for all s > sq = where a = a — b. If the nonlinearity F(c), which in 

this particular case equals c{a — c), satisfies the condition F'{c) < F'{0), then the minimal 
speed s is determined by F'{0) and not by the nonlinearity itself. In this sense the minimal 
speed is stable with respect to the perturbation of the nonlinearity. It is also known that 
the solutions of the parabolic equation with localized in space initial conditions converge to 
the wave with the minimal velocity [8]. 

Consider next the two-term approximation in ()8.ip . We obtain the equation 

| = ,(l_,e)g + c(.-.). (8.2) 

The corresponding stationary equation 

d{l — ^c)c + c((T — c) = 

has periodic in x solutions if 6cr > 1. Indeed, it can be written as the system of two first 
order equations: 

c' = p, p' = -$(c). 
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Figure 6: Stable equilibrium behind the wave 

where $(c) = c(cr — c)/{d(l — 7c)). Since = 0, $'(o") > 0, then there exists a family of 
limit cycles around the stationary point c = o", p = of this system. 

Thus, if 7 = there are travelling waves with the limits c = and c = a. Each of these 
values is a stationary solution of the evolution equation, the first one being unstable while 
the second one is stable. The minimal speed is determined by the linearized about c = 
equation. If 7 7^ 0, the same homogeneous stationary solutions exist but c = a can become 
unstable. Therefore the wave can connect the homogeneous solution c = at one infinity 
with the periodic solution at another infinity. However, it cannot be done for equation (j8.2j) 
because the coefficient of the second derivative will change sign. We need to consider the 
complete equation ()4.1|) . Its linearizations about c = and about c = a have the same 
properties as for equation ()8.2|) . Hence we can expect the existence of a periodic in time 
travelling wave, that is 

c(x, t) = u{x — st, t), 

where u{x,t) is a periodic in t solution of the corresponding equation. If we assume that 
the minimal speed is determined by the linearized about c = equation, as before, then it 
will have the same value sq = 2\fda. For the values d = 0.06, a = 2, 6 = 1 of the parameters 
considered in the simulations on Figure 6, we obtain sq ~ 0.48. This result is in a good 
agreement with the speed of propagation found numerically. 

It is interesting to compare the results on Figures 6 and 7. The values of the parameters 
are the same except for which determines the width of the support of the function 0. In 
the second case (Figure 7) it is twice more than in the first case (Figure 6). In agreement 
with the linear stability analysis the equilibrium behind the wave is stable in the first case 
and unstable in the second one. However the speeds of wave propagation in both cases is the 
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Figure 7: Periodic wave propagation 

same. This confirms the assumption that the minimal speed is determined by the system 
hnearized about c — 0, and that the solution converges to the wave with the minimal speed. 

9. Asymmetric evolution 

In Section 6 we considered the case where the function was even. Therefore its Fourier 

transform (p was a real-valued function, and the principal eigenvalue crossed the imaginary 
axis through the origin. On the other hand, the numerical simulations presented in Section 7 
showed that behavior of solutions in the case of the asymmetric function could be different 
in comparison with the symmetric case. Consider for certainty the following function: 

^^^^ \ , otherwise ^ ^ 

Then 

1 1 i 

m = ^ (cos(ez/) + isiYi{iy))dy = — sin(e7V) + — (cos(eA^) - 1). 

Hence the real part of the function $(^) is the same as before but the imaginary part is now 
different from zero. This means that the stability boundary is the same as for the symmetric 
function but the behavior of solutions is different. 

To describe the behavior of solutions at the stabihty boundary consider the hnearized 
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Figure 8: Asymmetric evolution 

equation 



dc d'^c 

— = d — - a I y)c{y)dy. (9.2) 



dt dx^ 
We look for its solution in the form 

c{x,t) = cos(^(x — st)). 

Substituting it in ()9.2|1 and separating the real and imaginary parts, we obtain 

/oo 
4>{y) cos{^y)dy = 0, 
-oo 

/oo 
sm{^y)dy. 
-oo 

The first equality is satisfied at the stability boundary. The second equality gives the relation 
between the speed s and the frequency ^. For the function (j9.ip we have 

' = ^ (cos(eiV) - 1). (9.3) 

We verify this relation with the numerical simulations shown in Figure 8. The values 
of parameters here {d = 0.01, a = 1,A^ = 1.1) are such that the homogeneous solution 
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is unstable but they are not far from the stabihty boundary (A^ = 1). The speed s of 
propagation equals the minimal speed Sq = = 0.2. The value of ^ taken from the 

numerical results is approximately 2.85. The right-hand side in ()9.3|1 equals approximately 
0.22. Therefore this equality is approximately satisfied. We note that it is satisfied at the 
stability boundary but not outside it. For larger values of the frequency of oscillations is 
quite different (see Figure 5, where all the parameters are the same as in Figure 8, except 
that N = 2) though the speed of propagation remains the same. 

10. Discussion 

In this work we study population dynamics with nonlocal consumptions of resources. 
We consider the conventional logistic equation with diffusion and with a nonlinear term 
describing reproduction and mortality of individuals in a population. The difference in 
comparison with the classical in the population dynamics model is that the death function 
depends not only on the density of the individuals at the given spatial point but also on their 
density in some its neighborhood. This is related to the assumption that the individuals in 
the population consume the resources not only at the point where they are located but also 
in some area around this point. 

The nonlocal consumption of the resources results in appearance of spatial structures. 
The homogeneous equilibrium, which would be stable otherwise, becomes unstable. This is 
a new mechanism of pattern formations different with respect to Turing structures. As it is 
well known, Turing or dissipative structures are described by reaction-diffusion systems of 
at least two equations for two different concentrations. They are often associated with an 
activator and inhibitor in some chemical systems. This mechanism can be briefly described 
as a short range activation and a long range inhibition. Though this mechanism is often 
applied to describe biological pattern formation, and in some cases it gives patterns very 
close to those observed in reality, it remains basically phenomenological (see [9]). 

The model considered in this work can show the emergence of patterns in a single pop- 
ulation, without competition between different species. To our knowledge, it is a new effect 
for modelling in population dynamics. Its biological interpretation is that an intra-specific 
competition can result in a spatial localization of the individuals in the populations. We can 
suppose that formation of herds of herbivorous animals or formation of local communities in 
human populations, though much more complex, can be related to this mechanism. 

The possibility of the appearance of spatial structures depends on the function (f){y). In 
particular, the piece- wise constant function introduced in Section 6 can lead to the instability 
but not the error function (Section 5). Let us discuss a possible biological meaning of the 
function (jfi.lj) . Consider an individual with his average location at y = where he has 
his family life (reproduction). He consumes resources during the day and returns to the 
family during the night. Suppose that he needs 8 hours per day to consume the sufficient 
amount of resources and that he can spend 12 hours per day outside the nest. Then he has 
2 hours to go to his work (consumption) place and 2 hours to return back. If we neglect the 
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Figure 9: Evolution of species according to Darwin 



cost of the displacement and if we assume that the working places are equally distributed, 
then he will go with a constant probability to any place in the two hours distance from 
the nest but not further than that. This gives us the piece- wise constant distribution ()6.1|) . 
If the diffusion coefficient in this case is sufficiently small, then a nonhomogeneous spatial 
distribution will emerge. We should distinguish here a fast motion home-work, which is 
considered as instantaneous displacement, and a slow motion due to diffusion. The latter is 
in fact a displacement of the nest. For example, it can be a relatively slow motion to move 
the home closer to the work place. The characteristic time of this motion can be related 
to the development of the infrastructure and of other factors. If the diffusion coefficient is 
large, then the spatial structure does not appear and the homogeneous equilibrium remains 
stable. In human populations it can result in appearance of big equally populated areas (e.g. 
Chicago area in the USA). 

On the other hand, if the cost of the displacement is high, then will decrease when 
\y\ increases. If the influence of this factor is sufficiently strong, then the population will 
also remain equally distributed though for another reason: to stay close to the consumption 
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place. 

Another context where the model considered in this work can be applied is evolution of 
species (see Section 3). In this case the space variable does not correspond to the physical 
space anymore but to the morphological space. Consider some morphological feature, for 
example, the size of bird's beak and the resources related to this feature. It can be grains of 
different sizes that birds consume. This brings us to the same piece-wise constant function 
(f) as before. By diffusion in this case we understand a slow change of the beak size from 
one generation to another due to random mutations or due to some other factors. The 
intra-specific competition for resources can result in the splitting of the population into one 
or several sub-populations (see Figures 2, 3) that can be considered as new species. If the 
diffusion coefficient is sufficiently small, then the new species emerge quickly in time and 
they are not connected to the previous species (Figure 4). This can be an explanation of the 
fact that transient forms in the evolution of species are not known. 

It is interesting to note that Figure 3 looks very close to the schematic representation of 
evolution given by Darwin in his "Origin of species" (see [10], pages 116, 117). 

Let us discuss finally the simulations shown in Figure 5. The function is not symmetric 
in this case. The biological interpretation of this assumption is that there is a preferential 
direction in the consumption of the resources. It can be related for example to a dominant 
direction of wind in the case of insects or to some other environmental or social factors (e.g. 
preferential transfer of jobs to the East and of workers to the West). In the evolution context 
it can be some asymmetry of morphological features. For example, if an animal can easier 
bent down to look for food than to stretch up it can result in appearance of species of a 
higher height. 
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